spatial analysis
Pascale
11/8/2019
libraries
library(sf)## Linking to GEOS 3.6.1, GDAL 2.2.3, PROJ 4.9.3library(dplyr)##
## Attaching package: 'dplyr'## The following objects are masked from 'package:stats':
##
## filter, lag## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, unionlibrary(ggplot2)
library(leaflet)
library(scales)
library(ggmap)## Google's Terms of Service: https://cloud.google.com/maps-platform/terms/.## Please cite ggmap if you use it! See citation("ggmap") for details.library(scales)
library(ggmap)
library(leaflet)prepare region shapefile
ak_regions <- read_sf("data/shapefiles/ak_regions_simp.shp")
plot(ak_regions)
datum is shape of sphere and projection is math used to flatten it
st_crs(ak_regions) # check projection## Coordinate Reference System:
## EPSG: 4326
## proj4string: "+proj=longlat +datum=WGS84 +no_defs"class(ak_regions)## [1] "sf" "tbl_df" "tbl" "data.frame"# switch projection to get rid of problem with international dateline
ak_regions_3338 <- ak_regions %>%
st_transform(crs = 3338) # crs is coordinate reference system
plot(ak_regions_3338)
summary(ak_regions_3338)## region_id region mgmt_area geometry
## Min. : 1 Length:13 Min. :1 MULTIPOLYGON :13
## 1st Qu.: 4 Class :character 1st Qu.:2 epsg:3338 : 0
## Median : 7 Mode :character Median :3 +proj=aea ...: 0
## Mean : 7 Mean :3
## 3rd Qu.:10 3rd Qu.:4
## Max. :13 Max. :4ak_regions_3338 %>% select(region)## Simple feature collection with 13 features and 1 field
## geometry type: MULTIPOLYGON
## dimension: XY
## bbox: xmin: -2175328 ymin: 405653.9 xmax: 1579226 ymax: 2383770
## epsg (SRID): 3338
## proj4string: +proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## # A tibble: 13 x 2
## region geometry
## <chr> <MULTIPOLYGON [m]>
## 1 Aleutian Islands (((-1156666 420855.1, -1159837 417990.3, -1161898 4169~
## 2 Arctic (((571289.9 2143072, 569941.5 2142691, 569158.2 214214~
## 3 Bristol Bay (((-339688.6 973904.9, -339302 972297.3, -339229.2 971~
## 4 Chignik (((-114381.9 649966.8, -112866.8 652065.8, -108836.8 6~
## 5 Copper River (((561012 1148301, 559393.7 1148169, 557797.7 1148492,~
## 6 Kodiak (((115112.5 983293, 113051.3 982825.9, 110801.3 983211~
## 7 Kotzebue (((-678815.3 1819519, -677555.2 1820698, -675557.8 182~
## 8 Kuskokwim (((-1030125 1281198, -1029858 1282333, -1028980 128403~
## 9 Cook Inlet (((35214.98 1002457, 36660.3 1002038, 36953.11 1001186~
## 10 Norton Sound (((-848357 1636692, -846510 1635203, -840513.7 1632225~
## 11 Prince William ~ (((426007.1 1087250, 426562.5 1088591, 427711.6 108999~
## 12 Southeast (((1287777 744574.1, 1290183 745970.8, 1292940 746262.~
## 13 Yukon (((-375318 1473998, -373723.9 1473487, -373064.8 14739~prepare the population data
pop <- read.csv("data/shapefiles/alaska_population.csv", stringsAsFactors = FALSE)
head(pop)## year city lat lng population
## 1 2015 Adak 51.88000 -176.6581 122
## 2 2015 Akhiok 56.94556 -154.1703 84
## 3 2015 Akiachak 60.90944 -161.4314 562
## 4 2015 Akiak 60.91222 -161.2139 399
## 5 2015 Akutan 54.13556 -165.7731 899
## 6 2015 Alakanuk 62.68889 -164.6153 777assuming WGS84 without metadata
pop_4326 <- st_as_sf(pop,
coords = c("lng", "lat"),
crs = 4326,
remove = FALSE)
head(pop_4326)## Simple feature collection with 6 features and 5 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -176.6581 ymin: 51.88 xmax: -154.1703 ymax: 62.68889
## epsg (SRID): 4326
## proj4string: +proj=longlat +datum=WGS84 +no_defs
## year city lat lng population geometry
## 1 2015 Adak 51.88000 -176.6581 122 POINT (-176.6581 51.88)
## 2 2015 Akhiok 56.94556 -154.1703 84 POINT (-154.1703 56.94556)
## 3 2015 Akiachak 60.90944 -161.4314 562 POINT (-161.4314 60.90944)
## 4 2015 Akiak 60.91222 -161.2139 399 POINT (-161.2139 60.91222)
## 5 2015 Akutan 54.13556 -165.7731 899 POINT (-165.7731 54.13556)
## 6 2015 Alakanuk 62.68889 -164.6153 777 POINT (-164.6153 62.68889)transform to match shapfile
pop_3338 <- pop_4326 %>%
st_transform(crs = 3338)calculate population by region
pop_joined <- st_join(pop_3338, ak_regions_3338, join = st_within)
head(pop_joined)## Simple feature collection with 6 features and 8 fields
## geometry type: POINT
## dimension: XY
## bbox: xmin: -1537925 ymin: 472627.8 xmax: -10340.71 ymax: 1456223
## epsg (SRID): 3338
## proj4string: +proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs
## year city lat lng population region_id region
## 1 2015 Adak 51.88000 -176.6581 122 1 Aleutian Islands
## 2 2015 Akhiok 56.94556 -154.1703 84 6 Kodiak
## 3 2015 Akiachak 60.90944 -161.4314 562 8 Kuskokwim
## 4 2015 Akiak 60.91222 -161.2139 399 8 Kuskokwim
## 5 2015 Akutan 54.13556 -165.7731 899 1 Aleutian Islands
## 6 2015 Alakanuk 62.68889 -164.6153 777 13 Yukon
## mgmt_area geometry
## 1 3 POINT (-1537925 472627.8)
## 2 3 POINT (-10340.71 770998.4)
## 3 4 POINT (-400885.5 1236460)
## 4 4 POINT (-389165.7 1235475)
## 5 3 POINT (-766425.7 526057.8)
## 6 4 POINT (-539724.9 1456223)total population by region
pop_region <- pop_joined %>%
as.data.frame() %>% # need to remove geometry because points that make polygon creates repition
group_by(region) %>%
summarise(total_pop = sum(population))
head(pop_region)## # A tibble: 6 x 2
## region total_pop
## <chr> <int>
## 1 Aleutian Islands 8840
## 2 Arctic 8419
## 3 Bristol Bay 6947
## 4 Chignik 311
## 5 Cook Inlet 408254
## 6 Copper River 2294pop_region_3338 <- left_join(ak_regions_3338, pop_region)## Joining, by = "region"plot(pop_region_3338)
pop_mgmt <- pop_region_3338 %>%
group_by(mgmt_area) %>%
summarise(total_pop = sum(total_pop), do_union=FALSE)
plot(pop_mgmt["total_pop"])
make maps!
rivers_3338 <- read_sf("data/shapefiles/ak_rivers_simp.shp")
st_crs(rivers_3338)## Coordinate Reference System:
## No EPSG code
## proj4string: "+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +units=m +no_defs"ggplot()+
geom_sf(data = pop_region_3338, mapping = aes(fill=total_pop))+
geom_sf(data = pop_3338, mapping = aes(), size=0.5)+
geom_sf(data = rivers_3338, mapping = aes(size = StrOrder), color="grey")+
scale_size(range=c(0.01, 0.2), guide = FALSE)+
theme_bw()+
labs(fill="Total Population")+
scale_fill_continuous(low = "khaki", high="firebrick", labels=comma)
incorporate basemap
pop_3857 <- pop_3338 %>%
st_transform(crs = 3857)# Define a function to fix the bbox to be in EPSG:3857
# See https://github.com/dkahle/ggmap/issues/160#issuecomment-397055208
ggmap_bbox_to_3857 <- function(map) {
if (!inherits(map, "ggmap")) stop("map must be a ggmap object")
# Extract the bounding box (in lat/lon) from the ggmap to a numeric vector,
# and set the names to what sf::st_bbox expects:
map_bbox <- setNames(unlist(attr(map, "bb")),
c("ymin", "xmin", "ymax", "xmax"))
# Coonvert the bbox to an sf polygon, transform it to 3857,
# and convert back to a bbox (convoluted, but it works)
bbox_3857 <- st_bbox(st_transform(st_as_sfc(st_bbox(map_bbox, crs = 4326)), 3857))
# Overwrite the bbox of the ggmap object with the transformed coordinates
attr(map, "bb")$ll.lat <- bbox_3857["ymin"]
attr(map, "bb")$ll.lon <- bbox_3857["xmin"]
attr(map, "bb")$ur.lat <- bbox_3857["ymax"]
attr(map, "bb")$ur.lon <- bbox_3857["xmax"]
map
}bbox <- c(-170, 52, -130, 64)
ak_map <- get_stamenmap(bbox, zoom = 4)## Source : http://tile.stamen.com/terrain/4/0/4.png## Source : http://tile.stamen.com/terrain/4/1/4.png## Source : http://tile.stamen.com/terrain/4/2/4.png## Source : http://tile.stamen.com/terrain/4/0/5.png## Source : http://tile.stamen.com/terrain/4/1/5.png## Source : http://tile.stamen.com/terrain/4/2/5.pngak_map_3857 <- ggmap_bbox_to_3857(ak_map)ggmap(ak_map_3857)+
geom_sf(data = pop_3857, aes(color = population), inherit.aes = FALSE)+
scale_color_continuous(low = "khaki", high = "firebrick", labels = comma)## Coordinate system already present. Adding new coordinate system, which will replace the existing one.
leaflet
custom function for alaska
epsg3338 <- leaflet::leafletCRS(
crsClass = "L.Proj.CRS",
code = "EPSG:3338",
proj4def = "+proj=aea +lat_1=55 +lat_2=65 +lat_0=50 +lon_0=-154 +x_0=0 +y_0=0 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs",
resolutions = 2^(16:7))pop_region_4326 <- pop_region_3338 %>%
st_transform(crs = 4326)leaflet(options = leafletOptions(crs = epsg3338)) %>%
addPolygons(data = pop_region_4326,
fillColor = "gray",
weight = 1)example 1 color palette
pal <- colorNumeric(palette = "Reds", domain = pop_region_4326$total_pop)
m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
addPolygons(data = pop_region_4326,
fillColor = ~pal(total_pop),
weight = 1,
color = "black",
fillOpacity = 1,
label = ~region) %>%
addLegend(position = "bottomleft",
pal = pal,
values = range(pop_region_4326$total_pop),
title = "Total Population")
mexample 2 markers sized by population
pal <- colorNumeric(palette = "Reds", domain = pop_region_4326$total_pop)
m <- leaflet(options = leafletOptions(crs = epsg3338)) %>%
addPolygons(data = pop_region_4326,
fillColor = ~pal(total_pop),
weight = 1,
color = "black",
fillOpacity = 1) %>%
addCircleMarkers(data = pop_4326,
lat = ~lat,
lng = ~lng,
radius = ~log(population/500), # arbitrary scaling
fillColor = "gray",
fillOpacity = 1,
weight = 0.25,
color = "black",
label = ~paste0(pop_4326$city, ", population ", comma(pop_4326$population))) %>%
addLegend(position = "bottomleft",
pal = pal,
values = range(pop_region_4326$total_pop),
title = "Total Population")
m